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ABSTRACT 

We explore the physics of shock evohition and particle acceleration in non-relativistic coUisionless 
shocks using multidimensional hybrid (kinetic ions/fluid electrons) simulations. We analyze a wide 
range of physical parameters relevant to the acceleration of cosmic rays (CRs) in astrophysical non- 
relativistic shock scenarios, such as in supernova remnant (SNR) shocks. We explore the evolution 
of the shock structure and particle acceleration efficiency as a function of Alfvenic Mach number and 
magnetic field inclination angle 0. We show that there are fundamental differences between high and 
low Mach number shocks in terms of the electromagnetic turbulence generated in the pre-shock zone 
and downstream; dominant modes are resonant with the streaming CRs in the low Mach number 
regime, while both resonant and non-resonant modes are present for high Mach numbers. Energetic 
power law tails for ions in the downstream plasma can account for up to 15% of the incoming upstream 
flow energy, distributed over ^ 5% of the particles in a power law with slope — 2±0.2 in energy. Quasi- 
parallel shocks with 9 < 45° are good ion accelerators, while power-laws are greatly suppressed for 
quasi-perpendicular shocks, 9 > 45°. The efficiency of conversion of ffow energy into the energy of 
accelerated particles peaks at 6* = 15° to 30° and Ma = 6, and decreases for higher Mach numbers, 
down to ~ 2% for Ma = 31. Accelerated particles are produced by Diffusive Shock Acceleration 
(DSA) and by Shock Drift Acceleration (SDA) mechanisms, with the SDA contribution to the overall 
energy gain increasing with magnetic inclination. We also present a direct comparison between hybrid 
and fully kinetic particle-in-cell results at early times; the agreement between the two models justifies 
the use of hybrid simulations for longer-term shock evolution. In SNR shocks, particle acceleration will 
be significant for low Mach number quasi-parallel fiows {Ma < 30, 9 < 45). This finding underscores 
the need for effective magnetic amplification mechanism in SNR shocks. 
Subject headings: cosmic rays — shocks — supernova remnants 
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Non-relativistic astrophysical shocks have long been 
considered as probable sources of acceleration of galac 



tic cosmic rays (C Rs) (Chen & Armstrong 1975 Bell 



1978||Blandford fc Ostrikcr,1978) . I'he CR spectrum ex- 
tends over many decades in energy, following a power law 
f{E) ex E~°' distribution. Nonthermal emission from 
supernova remnant shocks and from relativistic jets in 
Gamma-ray bursts and Active Galactic Nuclei is also 
modeled as synchrotron or inverse Compton scattering 
from a power law population of electrons at these sources. 
While the shape of the power law spectrum can be ex- 
plained by modern shock acceleration theory, particle ac- 
celeration efficiency, i.e., the number and energy fraction 
contained in the accelerated particles, depends on non- 
linear shock physics and has to be constrained through 
kinetic simulations. Such simulations are the subject of 
this paper. 
The first-order Fermi acceleration mechanism at shocks 



(geU 1978; Bla ndford fc Ostriker|1978||Blandford fc Eich- 
p^||1987| ) is the strongest contender to explain the ob- 
served power law spectra. The mechanism works by hav- 
ing particles cross a shock front several times, with a 
positive energy gain on every crossing. Diffusive Shock 
Acceleration (DSA) theory describes how these particles 
scatter upstream and downstream from the shock due 
to electromagnetic (EM) turbulence, and predicts a uni- 
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versal power law slope a, which to zeroth order depends 
only on the shock compression ratio {nd/uu, the ratio of 
the downstream to upstream density). Both the particle 
acceleration efficiency and the acceleration rate depend 
on the properties of turbulent fields and on the momen- 
tum of the accelerated particle; this is accounted for in 
the diffusion coefficient in the DSA model. Since the dif- 
fusion coefficient is not directly measurable in shocks of 
astrophysical origin, simplifying assumptions are usually 
made when estimating the maximum particle energies 
and acceleration times for a given object. One common 
assumption is to consider Bohm diffusion, by setting the 
diffusion coefficient to be proportional to the Larmor ra- 
dius of the accelerated particle. This assumption then 
constrains the energy gain to be directly proportional to 
time, E(t) ex t, and also limits the maximum energy of 
a particle by setting the maximum Larmor radius r^^ to 
be of the same order of magnitude as the system size. 
Another relevant acceleration mechanism is Shock 



Drift A ccelerat ion (SDA, Che n fc Armstrong 1975 Webb 
eral| [l983 Be gelman fc^r k 1990) and Shock Surfing 
Acceleration (S5X7Lee et al.||1996, |Hoshino fc Shimada 



2002 



Shapiro fc Uger 2003 ) , where a particle is acceler- 



ated by the E^ = —VxB upstream electric field in a 
magnetized shock, with B the upstream magnetic field, 
and V the upstream plasma fiow velocity measured in 
the downstream frame. The main difference between 
SDA and SSA is the importance of the electrostatic cross- 
shock potential for the confinement of particles in the ac- 
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celeration region; nevertheless, particles are accelerated 
by the same upstream electric field Eq, and for the pur- 
poses of this work we do not distinguish between the two. 
After being reflected at the shock, a particle will drift 
along the shock front in the direction of this upstream 
electric field, gaining energy. For magnetized shocks, this 
mechanism works along with DSA to produce energetic 
particles. 

Astrophysical collionless shocks can be characterized 
by their Alfvenic Mach number, defined as Ma = 
Vsh/VA, where Vgh is the shock front velocity (in the up- 
stream frame) and Va = Bo/^/innom is the Alfven ve- 
locity (with Bq, Uq and m being the upstream magnetic 
field, ion number density and ion mass, respectively). 
The magnetic field inclination angle 9 between B and 
the shock normal also plays an important role. Gen- 
eration of electromagnetic turbulence in the pre-shock 
zone is crucial for the DSA mechanism, since parti- 
cles need to be scattered back to the shock front mul- 
tiple times for efficient acceleration. Since cold par- 
ticles will be constrained to move along the magnetic 
field lines, different shock configurations will develop for 
angles 9 < 45° (quasi-parallel shocks), and for angles 
9 > 45° (quasi-perpendicular shocks). The main differ- 
ence for the two cases is that for 9 > 45° there is a sharp 
decrease in the number of particles escaping from the 
shock. For quasi-parallel shocks an extended pre-shock 
region will develop, with waves and turbulence being gen- 
erated by reflected ions that stream in front of the shock; 
quasi-perpendicular shocks will exhibit very little self- 
generated turbulence upstream, and instead a shock-foot 
region with the size comparable to the ion Larmor radius 
r^i will form. The DSA mechanism will then play a cru- 
cial role in particle acceleration in quasi-parallel shocks 
(for a compreh ensive review of microphysics of collision- 
less shocks see Treumaiin||2009 1 . 

Several approaches are commonly used to model the in- 
trinsically nonlinear processes that dominate shock acccl- 
eration. Semi-analytic kinetic models JKirk fc Heavens 



1989 



1991 1 [Achterberg et aL] KM_ 
tground tor the description of 



J Ballard & Heavens 
e a theoretical bac 
the particle acceleration problem. Similarly, Monte Carlo 
test particle simulations with prescribed turbule nce mod- 
els have been performed by several authors (e.g., Bednarz 
k Ostrowski|[!998; EUison fc Double 2004} [Niemlec fc 
Ostrowski 2004 1 . These can be extended over very large 



periods oi time and thus have predictive behavior over 
the long term evolution of particle acceleration. They 
lack the self-consistent properties of first principle simu- 
lations, such as fully kinetic particle-in-cell (PIC) mod- 
els. In PIC models, the full set of Maxwell's equations 
is solved on a finite grid, where particles are also pushed 
and serve as source terms for the creation of EM fields 
(Birdsall fc Langdon||1991 ). Hybrid simulations follow 



similar method and consider kinetic ions and fluid elec- 
since electron time-scales do not 



trons (Lipatov 
need to be reso 



2002D ; 

ved, evolution of the physical system for 



longer times is possible. 

Previous work on multidimensional PIC shock sim- 
ulati ons has been conducted for relativistic pair plas- 



mas (Spitkovsky"2005; •2008b[ JSironi fc Spitkovsky"2009[ 
[Nishikawa et al. 2009 ) , and foF relativistic ion-clcctron 
plasmas ( Spitkovsky|2008a| Martins et al.|2009 Sironi & 



Spitkovsky 20111). In the non-relativistic regime, both 
PIC (I Am'ano'ir Hoshino] 120071 J^ Takabe1[2008' 



Amano fc Hosliino||i>OOi)l 



J 



iqueime fc S pitkovsky 



201 1| ) and hybriJsimulations exist ([Eennett fc Ellison 



Ellison et al. 1993; Ciacalone et ai. 



1994, Giacalone et al. 1997' Scholer et al 



1992; Giacalone 
2002: Giacalon"e 



Electron accel- 



2004| Sugiyama[ 2011), aniong others 

eration is found to be more efficient in (magnetized) 

quasi-perpendicular shocks, in the non-relativistic case, 

while ion acceleration is mainly observed in quasi-parallel 

shocks. 

In this work, we concentrate on the systematic analy- 
sis of non-relativistic magnetized ion-electron shocks for 
diverse Mach numbers and magnetic field inclination an- 
gles. We extend previous work by analyzing a wide range 
of parameters, relevant for both solar wind shocks and 
for SNR shocks, and follow the shock evolution and ion 
acceleration for very long times and extended regions of 
space in two dimensions. Such ex tensive surveys are not 
commonly found in the literature ( Giacalone et al. ( 1997 1 
presents a similar survey, done in ID, and tor a more lim- 
ited range of Mach numbers). Qualitative differences are 
found between high and low Mach number shocks, and 
for quasi-parallel and quasi-perpendicular shocks. The 
analysis of individual ion trajectories provides a direct 
measurement of how particles are accelerated, and allows 
us to discriminate between the DSA and SDA models, 
and their relative contributions to the high-energy end of 
the particle distribution functions. We focus on particle 
acceleration and how the main shock parameters affect 
the acceleration efficiency. We also analyze the particle 
spectra behavior in time, and compare between spectra 
for the different Mach numbers and B field inclination 
angles, providing a complete survey of ion acceleration 
in non-relativistic magnetized shocks. 

This work is organized as follows. In Section [2] we 
present the simulation setup and the physical parame- 
ters to be explored. In Section [3] we show the global 
shock properties, including the main features of quasi- 
parallel and quasi-perpendicular shocks, together with 
typical particle spectra found in each case. In Section |4] 
we present the results of the complete shock parameter 
survey, analyze the shock properties with particular em- 
phasis on wave and turbulence generation, and its influ- 
ence on the measured particle spectra. In Section 5] we 
consider the particle acceleration mechanisms ancTdis- 
cuss how their relative contributions vary with the shock 
parameters. The flnal discussion, consequences and con- 
straints for astrophysical shock scenarios derived from 
the present work are presented in Section [6J Finally, we 
add an Appendix presenting a direct comparison between 
fully kinetic PIC simulations and their hybrid counter- 
parts, and we include a detailed analysis of the numerical 
parameters used. 

2. SHOCK SETUP AND SIMULATION PARAMETERS 

Two particle simulation codes are used to produce the 
results shown throughout this work: most results are ob- 
tained with the kinet ic ion/fluid electron code dHybrid 
(Gargate et al. 2007), and then checked wit h the rela- 



1993 

dimensional 



tivist i c tullylanencTiC code TRISTAN-MP (j Buneman 
"Spitkovskyjl2005|) 
tE 



m tins 
dimensional configuration 



While both codes are three 
work we only use them in two 



we only use tnem m 
The hybrid and PIC simu- 
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TABLE 1 
Key simulation parameters 



Run 


Ma 


e 


Box size (c/cjpi) 


A1...6 


3.1 


0,15,30,45,60,75 


2000 X 100 


B1...5 


6.0 


0,15,30,45,60 


2000 X 100 


C1...6 


10 


0,15,30,45,60,75 


5000 X 100 


V1...5 


31 


0,15,30,45,60 


8000 X 100 


£l..A 


100 


0,15,30,45 


8000 X 100 


T (PIC) 


4.7 





296 X 18 


Si 


4.7 





296 X 18 


Qi 


4.7 





2000 X 18 


Qz 


4.7 





2000 X 100 


Ga 


4.7 





2000 X 4 



lation models are very similar at the core. Both codes 
use a discrete grid where the electric and magnetic fields 
are calculated, which are used to advance particle ve- 
locities and positions in time; charges and currents from 
the particle populations are then used as source terms 
in Maxwell's equations to re-calculate fields. The hybrid 
approach differs from PIC in that the displacement cur- 
rent is neglected in Ampere's law, i.e.,V x B = AttJ/c, 
eliminating the propagation of light waves, and the elec- 
trons are modeled as a massless charge-neutralizing fluid. 
The electron equation of state is such that, after being 
adiabatically heated at the shock, the electron distribu- 
tion downstream will account for ~ 50% of the upstream 
ion flow energy; such a choice is physically motivated 
by results from PIC simulations. Details on the general 
properties of hybrid codes can be found iri Lipatov ( 2002 1 
and of PIC codes in Birdsall & Langdon (1991). 

In order to simulate a shock we consider an in-plane 
B field, propagate plasma through our simulation box 
in the —x direction, and refiect the flow from the left 
simulation wall. The two counter-propagating flows that 
result are electromagnetically unstable. The two flows 
thermalize and a shock propagating in the +x direc- 
tion is produced; the simulation frame is the downstream 
plasma frame for this setup. For electron-proton magne- 
tized shocks, which are the topic of this work, the initial 
plasma flow velocity V, the proton (or electron) number 
density n and the background magnetic field intensity 
and direction completely determine the physical proper- 
ties of the shock. For typical astrophysical shocks the 
inflowing upstream plasma is cold, such that the mag- 
netosonic Mach number is similar to the Alfvenic Mach 
number, and a strong (in the MHD sense) shock will form 
if Ma = V/Va > 1. 

The results are presented in normalized units, with the 
magnetic field normalized to the upstream background 
field Bq, the density normalized to the upstream plasma 
density no, the velocities normalized to the Alfven ve- 
locity, the spatial dimensions normalized to the ion in- 
ertial length c/ujpi (with ujpi — {Annoe'^ /m)^''^ the ion 
plasma frequency, and c the speed of light), and time 
normalized to the inverse ion cyclotron frequency u!~^ = 
{eBo/mc)~^. For hybrid runs, the grid cell size is set to 
0.2 c/u}pi, the time step is set to At = 0.01 oj^"'^, and 4 
particles per cell are used; total simulation times of up 
to 1200 a;^^^ were considered. For PIC runs, we use a 
grid cell resolution of 0.067 c/wp (electron skin-depths), 
and 4 particles per cell; a reduced rrii/me = 30 mass ra- 



tio is employed, in order to improve computation time 
efficiency. 

The full set of our runs is summarized in Table [ij 
Simulation box sizes vary with Mach number, from 
2000 X 100(c/wpi)2 up to 8000 x 100(c/wpi)2, for most 
of the hybrid runs. In the runs from sets A through £, a 
thorough physical parameter survey was conducted, in- 
cluding shock Mach numbers from 3.1 through 100, and 
B field inclination angles from 0° through 75°, as speci- 
fied in the table. In the following sections we thoroughly 
analyze the shock properties, particle spectra, and ac- 
celeration mechanisms observed, for runs in the sets A 
through £, and compare our results with previous simu- 
lations, in the frame of the DSA and SDA theories. Runs 
from the J- and Q sets were used in the direct comparison 
between PIC and hybrid results. Finally, extensive nu- 
merical convergence tests were also conducted for both 
codes. These tests included assessment of the effects of 
grid resolution, number of particles per cell, timestep, 
and numerical effects due to the finite box size (both in 
the longitudinal and transverse directions, as detailed in 
the Appendix). 

3. GLOBAL SHOCK FEATURES 

Figure [l] shows the transversely-averaged ion number 
density, B^ magnetic field component, x — Vx phase- 
space, magnetic field energy (normalized to the incom- 
ing flow energy), and the average magnetic fleld incli- 
nation angle (i.e., (arccos [i?a;/|-S|])) for Mach number 
Ma = 3.1 shocks for different B field inclination angles 
{e = 0°,30°,45° and 60°). The shock compression ratio 
is nd/riu ~ 4 in all cases, as expected from the MHD 
Rankine-Hugoniot relation^ From the phase space plot 
it is clear that for 9 = 60° there are essentially no parti- 
cles being reflected at the shock (Figure [It versus Figures 
fTfc,h and m). The properties of the magnetic field are 
very different: waves and turbulence are observable for 
quasi-parallel shocks {9 < 45°) in the pre-shock region 
upstream, and no magnetic field perturbations are seen 
for the quasi-perpendicular shock. These waves are gen- 
erated by the interaction of two counter-streaming beams 
in the pre-shock region: the upstream plasma and the 
particles refiected from the shock front. The waves are 
the source for the magnetic field turbulence that is neces- 
sary for the DSA mechanism to work. Their propagation 
direction is always parallel to the background magnetic 
field, such that the wave fronts reach the shock surface 
at an angle for all 9 > 0° (reading the propagation angle 
from Figure [T] is non-trivial, since the simulation box is 
not plotted with correct aspect ratio). 

As the cold upstream plasma reaches the shock, it 
heats up and starts to thermalize to a higher tempera- 
ture downstream. There is a temperature overshoot just 
after the shock jump, observable as a larger spread in 
the X — Vx plots in Figures [TJ:,h,m and r, where Tt > Td 
(with Tt and Td the temperatures in the transition and 
far downstream regions, respectively). The transition re- 
gion where Tt > Td also corresponds to a zone of en- 
hanced electromagnetic waves, which is observable as a 

^ This ratio depends on the adiabatic index which for 2D mag- 
netized case is the same as in 3D (particles' velocities are not con- 
fined to the simulation plane, and have 3 degrees of freedom in 
both cases). 
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Fig. 1. — Global shock features for 4 shocks, at time t = SOOui^^ , with Ma = 3.1 and = 0° (run yli, upper left 5 panels), 30° (run ^3, 
upper right 5 panels), 45° (run yi^, lower left 5 panels), and 60° (run C5, lower right 5 panels). For each run, from top to bottom we show 
the plasma density averaged over the y direction, Bz magnetic field component, x — v^ phase space, magnetic field energy normalized to 
the incoming flow energy (total in black, self-generated in red), and average magnetic field angle (i.e., the angle of the B field vector with 
the X direction in each cell, 8 = arccos [iJi/|_B|], averaged over the y direction). 



sharp rise in the magnetic field energy in Figures fTJi,i,n 
and s. FinaUy, Figures[T^,j,o and t also show the effective 
averaged B field angle along the shock direction. A grad- 
ual increase in effective 6 is most evident for the = 0° 
run. This feature of quasi-parallel shocks is due to the 
waves that are generated in the upstream pre-shock re- 
gion. The region of enhanced B field extends over a very 
large zone upstream in the 0° case, and it significantly 
affects the incoming plasma. Although the far upstream 
B field angle varies up to 60° for these runs, the mean 



B field angle at the shock boundary is nearly constant 
(~ 70°). This then suggests that even strictly parallel 
shocks are dominated by perpendicular shock physics, 
near the shock front, due to the strong transverse field 
of self-generated circularly polarized upstream waves. 

Figure [2] shows the global shock structure for Ma = 31 
shocks (the frames are the same as for Figure II]). In 
this case the shock compression is still ~ 4 as oefore, 
and the same suppression of refiected particles for an- 
gles 9 > 45° is observed. The most striking difference is 
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4000 
x|c/ai,i 



Fig. 2. — Same as FigurelT] for shocks with Ma = 31, and for the same inclination angles {d = 0° , d = 30°, 9 = 45° and 9 = 60°), at 
time t = 360 w^ . Waves in the upstream are affected by the limited transverse box size (see text). 



in the self-generated magnetic field: the waves excited in 
the upstream reach higher intensities relative to the back- 
ground magnetic field (this causes shock reformation in 
some cases, resulting in strong density fluctuations such 
as in Figurel2J-j at x ^ lOOOc/ojpi). Also, the shock tran- 
sition region is typically longer than for low Mach number 
shocks (- lOOc/wpi = 30rLi for Ma = 3.1 in Figure [l] 
compared to ~ 1500 c/wpi = 45 rLi for Ma — 31 in Figure 
[2]). The magnetic field energy profile for Ma = 3.1 and 
Ma = 31 shocks (red lines in Figure 111 and Figures l2li,i,n 
and s) is also quite different. For thenigher Mach num- 
ber shocks, the self-generated magnetic field energy (mi- 
nus the component due to compression, red lines) peaks 



at ~ 10% to 18% of the incoming plasma flow energy, 
whereas this figure is ~ 50% for the low Mach number 
case. Notice also that, although the Ma = 31 shock sim- 
ulations show the magnetic field propagating parallel to 
X, this is an artifact of the limited size in the y direc- 
tion; we used larger boxes to confirm that waves prop- 
agate parallel to the background B field, and that no 
other significant differences are observed in the results 
(see Appendix). 

Figures [3] and |4| show the detailed features for the 
quasi-parallel shock with Ma = 3.1 and 6 = 0°, and 
for the quasi-perpendicular shock with the same Mach 
number and 9 = 60°. In the quasi-parallel case, panels 
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Fig. 3. — Shock features for the Mj\ = 3.1 parallel shock {d = 
0°), run y4i, at i = 600 o;^ . Panel a) shows the plasma density, 
panel b) shows the density averaged over the y direction, panel c) 
shows the B^ magnetic field component, panel d) shows the lineout 
of the By (blue) and B^ (red) magnetic field components (along 
the dashed line), and e) shows the x — v^ phase space. Panel f) 
shows the particle spectra ~ 195 c/tjpi behind the shock (fitted to 
a Maxwellian, a power law, and an exponential cutoff). 

c) and d) show a clear signature of right-hand circularly 
polarized Alfven waves. In the simulation / downstream 
frame, the waves in the upstream are propagating toward 
the shock, parallel to the background magnetic field and 
with the B field vector rotating clockwise, when look- 
ing into the —x direction (wave is resonant with the ions 
traveling in the -\-x direction and k = — fc^. etj, see Figure 
Isli). The effect of these waves on the plasma is visible 
as density perturbations in the upstream region (Fig- 
ures pb,b). The particle spectrum shown in Figure pf 
(measured downstream at a; = 450c/a;pi) is fitted to a 
Maxwellian component (in red), and a power-law com- 
ponent with an exponential cutoff (in green). 



f{E) (X E^/h 



kiE- 



(1) 



where the power-law term is only added above an injec- 
tion energy Einj . The thermal spread for the Maxwellian 
distribution, expressed here and throughout the paper 
as a fraction of the incoming upstream plasma bulk-flow 
energy, is Eth = 0.39 Eup (39% of the upstream bulk en- 
ergy). From Rankinc-Hugoniot conditions the total en- 
ergy spread in the ion distribution should be ^ 0.5 Eup; 
from simulation results we instead find this value to be 
^ 0.42 Eup (which includes the 0.39 Eup in the thermal 
component), with the remaining energy (0.58 Eup) being 
stored mostly in the electron component. The power-law 
index measured in this case is a ~ 1.8 and is not changing 
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Fig. 4. — General shock features for the M^ = 3.1 quasi- 
perpendicular shock (9 = 60°), run ^5, at t = 600 w^ . Panel 
a) shows the plasma density, panel b) shows the density averaged 
over the y direction, and panel c) shows the x — Vx phase space. 
Panel d) shows the By (blue color table) and B^ (red color table) 
magnetic field components in the zoomed-in box (as depicted in 
panel a). Panel e) shows the lineout of the two field components 
over the y direction (along the green dashed line in panel d). Panel 
f ) shows the particle spectra ~ 195 c/ojp; behind the shock (fitted 
to a Maxwellian, a power law, and an exponential cutoff). 

significantly for t > 200 cj^^. There are, however, both 
spatial and temporal fluctuations in the spectral slope of 
up to ±0.2, if measurements are taken at different times 
and / or different regions downstream. The general trend 
shows a slightly steeper spectral index a ~ 2.2 at earlier 
times, when the shock is forming. These variations are 
mostly due to numerical effecttrl the energy efficiency for 
a given shock, however, does not vary significantly when 
numerical convergence studies are performed. 

The strong contrast with quasi-perpendicular shocks is 
evident in Figure El The upstream plasma density (panel 
a) is not perturbed and EM turbulence upstream is non- 
existent (panel d); also there are no particles returning 
upstream (traveling in the -l-a: direction). Electromag- 
netic waves are seen only in the shock transition region, 
and are propagating parallel to the background B field in 
this zone (Figures l4a,e). The resulting particle spectra 
can be seen in FigureK. Again, there is a thermal spread 
Eth = 0.39 Eup in this case, which is the value expected 
from the shock jump conditions for this B field inclina- 
tion angle. The non-thermal component of the spectrum, 
however, is very limited in range and the maximum par- 
ticle energy growth in time is negligible; the green line in 
the plot is still a fit to a power-law with index a — 2, but 

* Varying numerical factors such as resolution and number of 
particles per cell alters the noise properties of hybrid codes, which 
in turn influences the thermalization of the ions downstream, and 
injection into the power-law. 
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with a very steep cutoff starting at Ecut = 2Eup, which 
indicates that indeed particle acceleration is almost non- 
existent in this case. Similar results are found for run Ae 
for the same Mach number and 9 = 75°. 

4. SHOCK PARAMETER SURVEY 

The extensive parametric study conducted here in- 
cludes shock Mach numbers from Ma = 3.1 up to Ma — 
100, and both quasi-parallel and quasi-perpendicular 
shocks with angles from 0° up to 75°, as specified in 
Table IT] (runs A through £). The goal of this study is 
twofold: to understand how shock acceleration efficiency 
changes with these parameters and to understand the mi- 
crophysics of particle acceleration for the most relevant 
scenarios. In order to understand the particle accelera- 
tion mechanisms, it is also crucial to understand how EM 
turbulence is generated, which relates closely to particle 
diffusion and acceleration in the DSA framework. 

4.1. Wave generation and shock reformation 

The dominant waves present in magnetized non- 
relativistic shocks are electromagnetic ion modes. There 
are two relevant modes to consider: the right-hand 
Alfven mode, and the left-hand ion-whistler mode. The 
phase velocities for these modes in the hybrid limit are 



V, 



RH 



VaJi 






Vlh^VaJi 






(2) 



(3) 



These can be derived from the full wave dispersion re- 
lation for a cold magnetized ion electron plasma, in the 
non-relativistic limit of Va <Si c, considering a large ion 
to electron mass ratio, and typical frequencies w < uJci- 
Both full PIC and hybrid simulations show similar wave 
structures, and exhibit the modes described by Eq. ^ 
the agreement between the models is also very good 
see Appendix). 

These wave modes are seen in Figure |3] and Figure l4j 
As observed from the downstream simulation frame, the 
waves in the upstream are propagating towards the shock 
front along the background magnetic field; for the quasi- 
parallel Ma = 3.1 shock case (Figure Isl) the waves are 
right-hand circularly polarized Alfven waves with Vph = 

Vrh from Eq. m (i.e., when looking into the k direction 
of propagation at a fixed position in space, the B and E 
vectors turn clockwise in time). 

Figure [5] shows the averaged density and the lineout of 
the By and Bz magnetic field components at two different 
times for the same 9 — 0°, Ma = 3.1 run. It is apparent 
that the shock structure is identical at t — 400 w^"'^ and 
at i = 600 w^^, and in fact this is true also for larger 
times up to the end of the simulation. Although the 
pre-shock region continues to extend farther into the up- 
stream, the shock front propagates at a constant veloc- 
ity and the shock evolution is continuous in time. The 
zoomed-in boxes in Figure [5] also show that the largest 
wavelength in this upstream region (measured at a con- 
stant distance from the shock front) is A ~ 30c/a;pi, and 
it is similar for the two times considered. These waves are 
resonant with the highest-energy particles (ions) propa- 
gating away from the shock (traveling in the -{-x direc- 
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Fig. 5. — Time evolution for the averaged density (black), By 
(blue) and Bz (red) magnetic field components at times t = 
400a;^^ (panel a), and t = 600uJ~^ (panel b) (Ma = 3.1, = 0°, 
run Ai). Panels c) and d) show the magnetic field components at 
a fixed distance in front of the shock at the two times. 
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Fig. 6. — Time evolution for the averaged density (black). By 
(blue) and Bz (red) magnetic field components at times 4 = 60aj~ 
(panel a), t = 120a;^ (panel b), and t = ISOo;^ (panel c){Ma = 
31, 9 = 0°, run Vi). Panels d), e) and f) show the magnetic field 
components at a fixed distance from the shock front for the three 
time values. 
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tion); this can be inferred from both Figure^ and Figure 
[Sj by looking at the two components of tne magnetic 
leld. 

This behavior changes when we consider higher Mach 
number shocks. One to one comparison of the averaged 
density and magnetic field components for a, 9 — 0°, 
Ma = 31 shock can be seen in Figure [6] for three dif- 
ferent instants in time. In this case, the waves in the 
upstream zoomed-in boxes are still propagating towards 
the shock, but instead they are left-hand circularly po- 
larized ion- whistler waves with V^/j = Vlh from Eq. (3|), 
and they are non-resonant with the particles traveling 
away from the shock in the +x direction. A similar be- 
havior is seen for quasi-parallel shocks with Ma = 10 and 
Ma = 100. Also, shock evolution at later times in the 
high Mach number regime exhibits both the left-hand 
and the right-hand modes; which of the modes domi- 
nates depends on the local properties of the interacting 
flows, which vary with distance to the shock front - the 
resonant right-handed mode dominates on larger scales 
A ^ lOOOc/wpi, while the non-resonant left-handed mode 
has a typical wavelength A ^ lOOc/wpi. 

One very important consequence of the different be- 
havior between high and low Mach number shocks is vis- 
ible in Figure |6] for Ma = 31 shock. The self-generated 
magnetic field amplitude is greater than the background 
magnetic field by SB/Bq ~ 2 to 4 in some zones up- 
stream. This results in an effective magnetic field inclina- 
tion angle 9 > 45° in these zones, which in turn strongly 
affects the plasma flow. Figures [6^-c show strong den- 
sity perturbations growing in the upstream pre-shock re- 
gion, which eventually causes particles to reflect at these 
points and the shock front to reform. This discontin- 
uous behavior is in strong contrast with the results for 
shocks with lower Mach number; a quasi-steady behavior 
of the shock structure in the case of high Mach numbers 
is not guaranteed, although it might be achieved at later 
times. Another indication that quasi-stationary behav- 
ior was not yet reached at t — 180 w^^ is that the typical 
wavelengths observed in Figure |6}i-f are growing in time 
from A ^ 30 c/wpi at i = 60 w^ up to A ^ 75 c/wpj at 
t = 180 w^^. Notice, however, that the strong reforma- 
tion process in Figure |6] is a transient effect of the shock 
formation process, due to the initial low temperature of 
the two counter-propagating flows. Shock reformation 
later in time still occurs, but is more localized in space 
(shock reforms over distances ^ lOOc/ojpj). 

The strong non-resonant SB/Bq magnetic field gener- 
ation for high Mach numbers {Ma > 10) resembles the 
Cosmic Ray Current Driven (CRCD) instability. This 
electro magnetic instability was first proposed by |Bell| 
( 2004 1 in the context of magnetic field amplification in 



supernova remnant shocks, and further explored numer- 
ically in the non-linear phase by sev eral auth ors ("Ohiral 
et al. 2009 Stroman et al. 2009; Riq uelme fc Sp it kovsky, 
20D9[ [Gargate et al.||2010p .~ Within the framework of 
this mechanism, the CRs non-resonantly excite an elec- 
tromagnetic mode in the precursor of a SNR shock; the 
wavelengths are much smaller than the Larmor radius of 
the CR particles, and the polarization of the mode is the 
same as the ion-whistler wave branch described by Eq. 
([3|. In our simulations, however, the typical wavelength 
measured in the precursor (which is growing in time as 



the shock evolves) , is of the same order of magnitude as 
the Larmor radius of the lowest energy CR particles, and 
the waves propagate with a phase velocity close to Va- 

4.2. Shock acceleration efficiency 

We define shock acceleration efhciency in two ways: as 
the fraction of energy in the tail of the particle distri- 
bution function to the total energy in the distribution 
(energy efficiency), and as the number of particles in the 
tail of the distribution to the total number of particles 
in the distribution (number efficiency). Spectra in the 
simulations can be fitted to Eq. (IT]), and the injection 
energy Einj is then used as a threshold to define which 
particles are part of the tail of the distribution, or part of 
the thermalized Maxwellian component. To assure con- 
sistency of measurements between different shocks, we 
measure the spectra by considering a shell of particles in 
the downstream region just beyond the shock transition 
region (shell width is '^ 2 to 5rLi, and extends up to 
^ lOOc/wpi from the shock front for Ma — 3.1 shocks, 
and up to ~ ISOOc/wpi for AIa = 31 shocks). 

The particle spectra are also measured at the same 
time t — 300 o;^^ across all shocks. This time was chosen 
due to physical and numerical considerations. On the 
one hand, before t ~ 40 w^^^ the shocks are still forming, 
and thus measuring particle spectra at this time is not 
interesting. On the other hand, due to the finite size of 
the simulation boxes in the x direction, the production 
of energetic particles at the shock saturates after t ~ 
600 w^ for our numerical parameters; this is a known 
effect due to the fact that the most energetic particles 
start to escape from the simulation box, which introduces 
a high-energy cutoff in the spectrum. 

Figure [7] compares particle spectra measured for par- 
allel 9 = 0° shocks with different Mach numbers, for 
Ma = 3.1 (frame a). Ma = 10 (frame b), and Ma = 31 
(frame c), at time t = 300(u;^^. As time evolves, the 
maximum power-law tail energy grows and the temper- 
ature of the distribution decreases (Figure 13 1 in the 
Appendix shows this trend for a Ma = 4.7 shock). The 
energy efficiencies derived from fits as the ones in Figure 
[7] show a general decrease at larger Mach number; this is 
because there are fewer particles in the tail of the distri- 
bution (above the injection energy) as the Mach number 
increases. 

The shock acceleration efficiency as a function of mag- 
netic field inclination angle for the runs in sets A through 
V (Mach numbers from 3.1 to 31) can be seen in Figurejsj 
As a general trend, the acceleration efficiency decreases 
with increasing Mach number and with increasing angle. 
The highest energy efficiency measured was E^ff ^ 15% 
for the Ma = 6, 9 = 15° shock. It is also apparent from 
Figure J8] that there is a consistent increase in efficiency, 
for eacn Mach number, at an angle of 15° to 30°. There 
is thus a critical intermediate angle for which the shock 
energy conversion efficiency is greatest. This can be un- 
derstood in the framework of DSA and SDA theories: 
the SDA mechanism, which is faster than DSA, is more 
effective for higher B field inclination angles. 

The energy efficiency decrease with shock Mach num- 
ber is very clear from Figure [8] The trend is similar when 
looking at number-efficiency (fraction of particles in the 
tail to the total number of particles), which is seen in 
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Fig. 7. — Particle spectra (and function fits) measured in the 
downstream at time t = 300 ix)^^ (black crosses). Spectra are 
shown for run A\ (M^ = 3.1, 9 = 0°) in frame a), for run Ci 
{Ma = 10, e = 0°) in frame b), and for run Di {Ma = 31, 6» = 0°) 
in frame c). 
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Fig. 8. — Acceleration efficiency, defined as the fraction of energy 
in the power law tail to the total energy in the distribution, as a 
function of the shock magnetic field inclination angle. Efficiencies 
for runs Ai to Ai {Ma = 3.1) are in blue, runs Bi to B4 {Ma = 6) 
in red, runs Ci to C4 {Ma = 10) in green, and runs ©i to 1)4 
{Ma = 31) in grey 

Figure [71 where the power law index is roughly the same 
a ^ 2±0.2 for all runs at t = 300uj~^^, but the number of 
particles in the tail decreases from frame a) {Ma = 3.1) 
to frame c) (M^ = 31). This result is an indication that 
particle injection into the DSA process is less efficient at 
higher Mach numbers due to the local conditions near 
the shock front. Analyzing the micro-structure of the 
shock transition, and understanding in detail the parti- 



cle acceleration mechanism is thus of crucial importance. 
It should also be noted that there is a a slight overall in- 
crease in the efficiency from Ma = 3.1 to Ma = 6 shocks 
(Figure Is] blue line to the red line) which is c ontrary to 
the general t rend; a similar increase is found in [Giacalone| 
et al.| ( |1997|) at these Mach numbers. This is due to two 
effects : on the one hand there is always an increase of the 
maximum energy in the ion distribution function as the 
Mach number increases, but on the other hand there is a 
consistent decrease of the number of accelerated particles 
in the tail for higher Mach numbers. Since a '^ 2, there is 
an equal amount of energy per decade in the distribution 
and thus high energy particles contribute significantly to 
the energy efficiency. The two effects are then both rele- 
vant. At low Mach numbers from Ma = 3.1 to Ma = 6, 
the first effect dominates, whereas for Ma = 10 and be- 
yond the second effect dominates, which then results in 
the overall efficiency decrease seen in Figure |8] 

For higher Mach number shocks {Ma — 100), the 
general shock features, such as the electromagnetic field 
structure, including the polarization of the waves in the 
upstream and the shock compression ratio are consis- 
tent with the shock jump conditions (and are similar 
to Ma = 31 shocks). However, in this case, we found 
that shock evolution would have to be followed for longer 
times and using a longer simulation box for a consistent 
measurement of the efficiency; tentative values in the or- 
der of ~ 1% were found at the end of the simulations 
at i = 80 w^^ for A4a = 100 shocks. For shocks with 
angles 6 > 45° we found that the power law was strongly 
suppressed. 

5. PARTICLE ACCELERATION 

As the particle acceleration and the shock evolution 
are dynamic processes, it is expected that the maximum 
particle energy in the power-law tail should grow in time. 
For a DSA-accelerated particle, the total elapsed time 
for a particle to accelerate from momentum po to p de- 
pe nds on the diffusion coefficient and can be estimated 
as prury|ri983) 



M 




(4) 



where Vu, Vd are the upstream and downstream fiow ve- 
locities in the shock frame, and k„, Kd are the upstream 
and downstream spatial diffusion coefficients along the 
shock propagation direction. If a Bohm-type diffusion 
coefficient is now considered, n = 1/3 Aw, where A = tli 
is the particle's mean free path and v the particle's ve- 
locity, it is possible to deduce Tacc oc E; that is, under 
Bohm diffusion, the energy of a particle undergoing dif- 
fusive shock acceleration will grow linearly with time. 

Figure [9] shows the evolution of the energy of the most 
energetic particle in a given simulation for shocks with 
different B field inclination angles. Energy evolution 
in time is fastest for the Ma — 3.1, 45° shock, which 
also corresponds to the case where the largest energy 
is attained. This is partly due to the SDA mechanism: 
the upstream background electric field intensity will be 
higher for higher inclination angles. As it can be seen in 
Figure Is] (blue line), however, the Ma = 3.1, 45° shock 
is not tne one with the best acceleration efficiency; this 
is due to the fact that injection into DSA mechanism 
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Fig. 9. — Evolution of maxiinuirL particle energy in time, for 
shocks with Mj^ = 3.1 (runs from set yi), for angles 6 = 0° (solid 
black line), 9 = 15° (red), 9 = 30° (green), 9 = 45° (blue), 9 = 60° 
(dashed line), and 9 = 75° (dotted line). 
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Fig. 10. — Time evolution of key characteristics for a typical ac- 
celerated particle, extracted from run ^i (M^ = 3.1, 9 = 0°). 
From left to right wc show: a) particle energy, b) particle distance 
from the shock (in the shock propagation direction), c) the product 
of the Energy and the particle's velocity components E Vx (red), 
Evy (green), and Ev^ (blue), and d) the displacement of the parti- 
cle in the x (red), in the y (green) and the z (blue) direction. Panel 
b) also shows the B^ magnetic field component time evolution at 
the shock front as color. 

is harder for higher inchnation angles, such that a lower 
number of particles are present in the power-law tail. The 
evolution of energy in time for the Ma = 3.1 shocks with 
intermediate inclination angles 6 ^ 15° to 30° is slower 
than for the same shock with 9 — 45°, and faster than 
for lower angles. Since injection into the DSA mechanism 
is more effective at lower angles, the end result is that 
shocks with angles in the range of ~ 15° to 30° are the 
most efficient particle accelerators, for the parameters 
surveyed. 

By measuring the maximum energy in the particle dis- 
tribution at several points in time in one simulation, it 
is possible to estimate what is the energy gain in time. 
For the runs with AIa = 3.1 we find that the energy evo- 
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Fig. 11. — Time evolution of key characteristics for a typical 
accelerated particle, extracted from run A.4 {M_a = 3.1, 9 = 45°). 
Quantities plotted are the same as in Figure [To] 

lution follows a power- law in time, a,s E oc t^ , except for 
angles 9 > 60, which show almost no energy growth in 
time. For the Ma = 3.1, e = 0° run, we get l3 ~ 0.28, 
which implies that diffusion is less efficient than Bohm 
(i.e., using Eq. (|4]), /3 < /3b = 1), in the approximation 
that all the energy gain is due to DSA. This measurement 
yields similar results for other B field inclination angles; 
a systematic study of the variation of /?, and of the diffu- 
sion coefficient itself, with Mach number and inclination 
angle will be the subject of future work, since this study 
would imply a different and more thorough analysis. 

The trajectories of energetic simulation particles sup- 
port the claim that the DSA mechanism is responsible 
for particle acceleration in non-relativistic magnetized 
shocks. Particle track information is presented in Figures 
10 and 11 for typical high-energy particles in Ma = 3.1 
sEocks with 9 = 0° and 9 = 45°, respectively. These 
resemble th e particle t r acks a nd acceleration mechanism 
described in Sugiyama ( 2011[) , where particles are greatly 
accelerated only within an effective distance ^x^ff from 
the shock front. 



Figure 10 1 shows the energy increase with time, from 
around the injection time at i ^ 1000 cj^ , which corre- 
sponds to the first crossing of the shock, until the par- 
ticle escapes into the upstream medium. Multiple up- 
stream/downstream reflections can be observed in Fig- 
ure [TOJa, confirming a diffusive- type acceleration process. 
Most of the energy gain for this representative particle is 
achieved within a distance Axe// ~ 30 c/wpi of the shock 
front, where the electromagnetic waves are more intense, 
and thus are more effective in scattering particles back to 
the shock. Also, this energy gain was achieved mostly in 
the time interval from t ^ 1090uj~^^ up to i ~ llSOw^'^, 
and was approximately linearly proportional to time. For 
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this particular shock case, where 6 = 0° far upstream, 
particle confinement near the shock front is greatly en- 
hanced by the self-generated EM waves which effectively 
increase the averaged magnetic field inclination angle in 
a region Ax^p ^ 100 c/copi upstream of the shock front 
(see Figure lie). Thus, although the energy gain for in- 
dividual particles is approximately linearly proportional 
to time over finite time intervals, the overall energy gain 
for the ensemble of accelerated particles over longer time 
intervals follows E oc i"'^^ because, most of the time, the 
particles are not in the Ax^ff upstream region where 
acceleration is effective. 

In the case of the Ma = 3.1, 9 = 45° shock. Figure 
11 accelerated particles also undergo several reflections 
upstream, and are carried back to the shock. In this case, 
the particles are effectively reflected back to the shock up 
to larger distances Ax^p ^ SOOc/wpj. Since the power 
and polarization of the EM waves generated due to the 
shock is similar between the 6 = 45° and 6* = 0° cases, the 
increase in the characteristic length for effective reflection 
is just due to the higher B field inclination angle (i.e., 
particles are more likely to be carried back to the shock 
when 9 is higher). 

Overall, the global physical picture of the dependence 
of ion acceleration on the magnetic field angle is as fol- 
lows. The upstream cold flow streaming towards the 
shock will scatter at the shock front, and there will be 
fewer particles escaping back upstream with increasing 
angle. This greater injection efficiency for lower angles is 
theoretically de scribed in the literature (e.g., iGiacalone) 
& Jokipii 1996): for a given velocity and B field inten- 
sity, a larger angle implies a larger B field component 
parallel to the shock front, which will prevent an in- 
creasing number of particles from streaming away from 
the shock in the upstream direction. Once a particle 
is injected into the acceleration mechanism, however, a 
larger 9 will keep higher energy particles closer to the 
shock front where they will be accelerated faster. Our 
results show just that: the typical acceleration rate for 
a 6* = 0° case is AE/At - O.eEupWci (Figure [lO| from 
t ~ 1090 to t - 1180, w^^); for a 6* = 45° shock this rate 
is AE/At - 5Eupa;ci (Figure [n] for t > 750^^1). The 
acceleration of individual particles is due to the upstream 
convective electric field (which exists even for = due 
to upstream waves), and thus the energy gain is directly 
proportional to the electric field strength and to time. 
The overall particle acceleration is slower {\.e.,E ex f^, 
with /3 ~ 0.28), and is only slightly faster for larger angles 
due to the increasing intensity of the convective electric 
field. 

Particle injection into the acceleration mechanism also 
changes with the shock Mach number. By analyzing a 
large number of trajectories, for both high and low Mach 
numbers at angles 9 < 45°, it is apparent that almost all 
particles are injected into the acceleration mechanism by 
reflection at the shock front, instead of thermal leakage 
from the downstream side. This thermal leakage pro- 
cess, which depends on the micro-physics of the shock 
structure, i s a parameter in Diffusion Shock Accelera- 
tion t heory( |Malkov fc Drury|2001[ |Amato fc Blasi||2005 
2006 1 , and is assumed not to vary with shock Mach num- 
ber. Our simulations, however, show that the number of 
particles that are actually injected into the DSA mech- 



anism decreases with Mach number, so that the number 
efficiency varies from Neff - 4.5% for Ma = 3.1, 6i = 0° 
down to Neff ~ 1% for Ma = 31, 9 = 0°, and even lower 
Neff - 0.12% for Ma = 100, 61 = 0° ahhough measured 
at an earlier time. 

The physical reason for this decrease in injection ef- 
ficiency with Mach number can be understood by an- 
alyzing the magnetic field energy at the shock front. 
In the shock transition region this energy, expressed 
as a fraction of the incoming fiow energy, will increase 
relative to the far upstream value, due to the waves 
present in this zone (see for example, the B field pro- 
file in the shock transition region in Figures [5] and [61 
or the 2D profile of the same waves in Figures [l] and 
^. However, the average fraction of energy in these 
waves, relative to the upstream flow energy, decreases 
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50% for the M . = 3.1, 6* = 0° case, 



down to AE^ave/Eup ~ 10% for the Ma = 31, 6* = 0° 
case, and even lower AE^ave/Eup ~ 2.5% for Ma = 100, 
9 = Q° . Also, for lower Mach numbers a sharper mag- 
netic barrier is observed at the shock, while for higher 
Mach numbers the lower AE^ave/ Eup magnetic barrier 
extends further downstream. Or, stated in a different 
way, the shock transition region is longer (measured in 
ion skin depths or as a fraction of the upstream Larmor 
radius) for higher Mach number shocks. Typical cold up- 
stream particles are then more easily specularly-reflected 
for low Mach numbers than for higher Mach numbers, 
which results in a larger fraction of particles being in- 
jected into the DSA mechanism for low Mach numbers. 

6. DISCUSSION AND CONCLUSIONS 

In this paper we have discussed ion acceleration in non- 
relativistic magnetized shocks. In an effort to understand 
particle acceleration, we conducted a complete 2D sur- 
vey of this class of shocks, varying both the Alfvenic 
Mach number and the magnetic fleld inclination angle 
for each shock. As shown before by other authors, our 
results confirm that there are fundamental differences 
between quasi-parallel and quasi-perpendicular shocks. 
Our study allows to systematically trace these differ- 
ences as functions of the fiow parameters and see the 
transition between different regimes. For quasi-parallel 
shocks, a significant number of particles are reflected by 
the shock back into the upstream medium, which causes 
electromagnetic instabilities that propagate back to the 
shock, and create a turbulent upstream medium. In con- 
trast, for quasi-perpendicular shocks, with 9 > 60°, there 
is no significant particle population being reflected at 
the shock front, and hence EM turbulence upstream is 
greatly reduced. 

Our results show a fundamental difference between 
shocks at low Ma < 10 Mach numbers, and shocks 
with Ma > 10; for low Mach number shocks a right- 
hand mode with SB/Bq < 1 is excited in the up- 
stream, resonant with the ions streaming away from the 
shock, whereas for higher Mach numbers a left-hand non- 
resonant mode predominates. In the high Mach number 
case, the left-hand non-resonant waves exhibit B field 
amplitudes SB /Bo > 1; these waves resemble the CRCD 
instability, as they have the same polarization. Here how- 
ever, the wavelength A is comparable to the Larmor ra- 
dius of the lowest-energy CR particles, and the phase ve- 
locity of the waves is around the Alfven velocity. Thus, 
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from the point of view of the highest-energy CR particles, 
since 5B/Bq > 1, X < r^i and Va <C Vjiow (i.e., the phase 
velocity of the waves is much less than the upstream flow 
velocity), the ion- whistler waves observed should have 
similar scattering properties to CRCD-generated electro- 
magnetic turbulence. 

For high Mach number shocks the self generated EM 
waves with 6B/Bq > 1 in the upstream medium start 
to reflect particles locally, and the shock is observed to 
reform further upstream; a similar effect of shock refor- 
mation due to Short Large-Amplitude Magnetic Struc- 
tures (SLAMS) has also be en described in the literature 
( Dubouloz fc Scholer|1995 ). This shock jump is observed 
at a relatively early stage t ^ 100 w^^ of the shock prop- 
agation, and it is thus possible that this is a transient 
effect due to the shock formation process. There are in- 
dications, from longer simulations (in time) over larger 
spatial domains that this is the case: although reforma- 
tion still happens at late times, its effects are less notice- 
able than at the time the shock is forming. These sim- 
ulations also show that the dominant EM modes change 
in time if the particle flux changes as well, which leads to 
a highly dynamic shock structure where sub-shocks form 
in localized regions of the upstream medium. Thorough 
analysis of such effects will be the subject of a future 
publication. 

We compared shock acceleration efficiencies in simi- 
lar conditions for all Mach numbers and angles, as de- 
scribed in Table [l] and showed that the highest energy- 
conversion efficiencies are achieved for intermediate an- 
gles 9 ^ 15° to 30° and low Mach numbers. The power 
law indexes are similar in most cases, yielding a = 2±0.2. 
The typical growth of maximum energy in time is found 
to be'E' ex t°-28 (for Ma = 3.1 and = 0°). This value is 
indicative of a diffusion process less effective than Bohm 
diffusion, which would correspond to E' ex i; a full com- 
parison of diffusion coefficients for different Mach num- 
bers and inclination angles will also be the subject of 
future work. 

The analysis of representative particle tracks for two 
quasi-parallel low Mach number shocks is consistent with 
DSA mechanism, and confirms that particles are re- 
flected multiple times towards the shock front. Over- 
all, we can conclude that ion acceleration in magnetized 
non-relativistic shocks is efficient at finite intermediate 
angles, in the quasi-parallel regime, and for low Mach 
numbers {Ma < 10 ). This is explainable if we take 
into account that injection into a diffusive acceleration 
process is more likely at low angles, since a significant 
number of particles are allowed to travel farther into the 
upstream medium, and generate EM waves and turbu- 
lence which enhance further particle reflection towards 
the shock. For higher inclination angles, however, a par- 
ticle with a sufficiently high energy will be more easily 
confined to the shock front, and will be able to accel- 
erate faster due to stronger shock drifting (Figure [9]); in 
terms of energy-conversion efficiency an intermediate an- 
gle yields the best results. Also, although higher Mach 
number shocks are able to accelerate particles to higher 
energies (as a fraction of incoming beam energy), the 
shock efficiency is lower because injection into the accel- 
eration mechanism is less efficient at higher Mach num- 
bers. 



It is interesting to compare the present results with 



a similar shock survey done in Giacalone et al. (19971 



(henceforth G97) using ID hybrid simulations tor typi 
cal solar wind flow parameters. For a more limited range 
of Mach numbers, from Ma = 1.5 up to Ma — 12 the au- 
thors also show downstream spectra of accelerated par- 
ticles; the power-law tails extend over a similar range 
in energy as in our case, although high-energy particle 
splitting, and the introduction of external EM turbu- 
lence was used to facilitate the diffusive process in the 
case of G97. In the comparable range of Mach numbers 
(up to Ma = 12), the conclusions about shock accelera- 
tion efficiency are similar. At higher Mach numbers (not 
explored in G97) we observe a clear decrease in efficiency, 
and different features such as shock reformation and the 
left-hand wave polarization mode. 

Energetic arguments, in order to explain CR energies 
of up to lO^'^ eV from SNR sources, usually require a large 
fraction > 30 — 40% of the total energy in the shock to 
be transferred to accelerated particles, which indicates 
that shocks with Ma < 30 and 9 < 30° are the most 
likely candidates for CR acceleration. For angles greater 
than 45°, injection is very limited and no significant ac- 
celeration is observed; it should be noted, however, that 
if pre-accelerated particles could be injected in shocks 
with 9 ~ 60° in sufficient numbers, it would be possible 
to have higher acceleration efficiencies due to enhanced 
shock drift acceleration in the shock front (although this 
does not seem to occur naturally in the simulations). 

Relatively low Mach numbers {Ma < 30) are thus re- 
quired for effective ion acceleration, which then implies 
that strong magnetic field amplification should occur at 
SNR shock sites (if we consider typical outflow velocities 
of ^ 1000 km/s in interstellar fields of 1 /iG). This strong 
B field am plification is inferred from observations (for a 
review see 



Reynolds et al 
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Garg ate et al.|2010p . Such low Mach numbers are also fa 
vorec TBy the electron accelera tion mechanisms in shocks 
([Riquelme & Spitkovsky 2011). We do see strong B field 
amplification in our current simulations, but the effective 
Mach number near the shock front is at most reduced by 
half (effective Mach number near the shock is Ma ~ 15, 
corresponding SB ~ 2Bo, for the case of Ma = 31, 6* = 
far upstream conditions). A magnetic field amplification 
of at least an order of magnitude would be necessary in 
a SNR shock. As such, a self-consistent shock simulation 
where both accelerated particles and the CRCD instabil- 
ity are observed and where SB ~ IOBq is still lacking. 

The numerical properties of the hybrid and PIC codes 
must be very well understood in order to validate results. 
We show the effects of dimensionality, and validate our 
choice of numerical parameters in the Appendix. As fu- 
ture work it would also be interesting to compare the 
current results with 3D simulations, in particular regard- 
ing the diffusion pro perties which can be different in 3D 



(Jokipii et al. 1993 



and the characteristics of the wave 
spectra generated for quasi-parallel shocks, both in the 
low and high Mach number regimes. Finally, a complete 
statistical study of self-consistently accelerated particles 
in hybrid and PIC simulations would be important in or- 
der to assess the dependence of the diffusion coefficient 
on shock properties (such as the Mach number and B 



Ion acceleration in shocks 



13 




□ 1 



_ 15 
/ 10 

E 



d) 



r(, 




-I). 



Ui 




■ft*^ 



m) 



N^ 



a 





50 


100 


150 

X [0/OJ,J 


200 


250 




50 


100 


150 
X [o/to J 


200 


250 




60 


100 


150 
x|c/«,„| 


200 


250 






























10° 

in' 








^ 


\1/ 


\- 


i) 




— 


A 


\1/ 


>» 


'tft 


0) 


^ 




7\: 




f<telffrl| ■ 


«) 


"^ 


"^ 




_y\ 


^ 


10* 
10* 

















E|1/2mv,,] 



E[1/2m v,„'| 



E [1/2 m v,h" 



Fig. 12. — Global shock features for a parallel {9 = 0°) shock with Mach number Ma = 4.7 at t ~ 87aj^j^ for a Tristan-MP fully kinetic 
run (J-') (left panels), an hybrid run (Qi) with the same numerical parameters (middle panels), and an hybrid run {G2) with an extended 
box size in the x direction. From top to bottom the frames show: the plasma density, the plasma density averaged over the transverse 
direction y, the Bz component of the magnetic field, Bz averaged over the y direction, and the particle spectra at three different positions 
{x = 40c/(.i;pi in black, x = 120c/i.i;pi in red, and x = 200c/ajpi in blue). 



field inclination angle), and also on the individual parti- 
cle energies. Such a study could then be used to identify 
the astrophysical relevant scenarios for the acceleration 
of cosmic rays. 
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APPENDIX 

COMPARISON OF PIC AND HYBRID RESULTS 

Comparing two different simulation models is useful to determine the domain of validity of t he models. This kind 
of analysis has been done in the context of shocks for ID hybrid and Monte Carlo simulations in Ellison et al. (1993) 



and for 2D hybrid and PIC simulations of strictly perpendicular shocks in [Hellinger et a~ ( 20U7p . Here we focus 
on comparing the hybrid and fully-kinetic methods for strictly parallel shocks, and in particular on how the shock 



on comparing the hybrid and fully-kinetic methods for strictly parallel shocks, and in particular on how the shock 
structure and measured particle spectra is affected by the models and the choice of numerical parameters. 

We chose to compare a Ma = 4.7, 6* = 0° shock. The agreement between the two models regarding the general 
features of the shock is very good, as it would be expected, for similar numerical parameters. The shock compression 
ratio is nd/n,j ^ 4, propagating at T^ ~ 1.6 Va in the downstream reference frame; also in agreement with the Rankine- 
Hugoniot conditions for a parallel magnetized shock. Figure [12] shows the ion density and B^ magnetic field component 
for runs J^, Qi and Q2 (from left to right, the top 4 panelsJTwith all quantities plotted using the same scales. For 
runs A and Bi the numerical parameters are the same (namely, the same box size 296c/wpi x ISc/wpj is used), and 
both the shock structure and the waves in the shock transition region are very similar. However, when we consider a 
larger L^ — 2000c/wpi box size for run Q2 (right panels in Fig 12), the size of the shock transition region is reduced 
significantly. The difference in the wave structure in the shock transition region, in this case, is due to the fact that the 
pre-shock region is not fully resolved for the smaller box size. This zone is strongly affected by the particles reflected 
at the shock front and traveling in the upstream medium towards +x; by limiting the box size to L^ = 296c/a;pi in 
T and Qi, the pre-shock zone is also limited. This limitation in turn affects the incoming particles (in the upstream, 
traveling towards the shock), which in the realistic case of run Q2 start to gyrate before reaching the shock, and thus 
forming a less coherent beam (which leads to less coherent waves). 

We show the measured particle spectrum at a; = 40 c/wpj in the downstream, x = 12 c/wpi in the transition region 

The small differences between 



12 



and X — 200 c/wp; in the upstream for the three runs in the bottom panels of Figure 
the spectra measured in the three runs A, Bi and B2 are due to the differing number of particles used in plotting the 
spectra (less particles are being used to plot the spectra in the hybrid case). The downstream spectra (black lines in 
the bottom frames of Figure IM can be modeled by a Maxwellian with a thermal spread of Etn ^ 0.3 Eup (^ 30% of 
the incoming flow bulk energy), a power law tail f{E) ex E~'^ and an exponential upper energy cutoff. 



The time evolution of the downstream particle spectrum can be seen in Figure 13 Panel a) corresponds to the hybrid 
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1.0 

Fig. 13.— Ion spectra for a Ma = 4.7, 9 = 0° shock, at times t = 40a;^^ (black), t = 87 wT^ (red), t = 200uj~^ (green), and t = 400 oj^^ 

(blue) for hybrid run Q-j (panel a). Panel b) shows the spectra at the same times (t = 40 w^ in black and t = 87(.j|rT in red) for the PIC 

run J^. Panel c) shows the hybrid spectra from run Q-i, at time t = 400 w^ fitted with a Maxwellian, a power law, and an exponential 
cutoff. 
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Fig. 14. — Global shock features for an hybrid parallel (9 = 0°) shock with Mach number M^ = 4.7 at t ~ uj^^ (run Q^). Numerical 
parameters are the same as in run Cy2, but an extended y dimension is used. From top to bottom, the panels show: a) plasma density, b) 
plasma density averaged over y, c) Bz component of the magnetic field, d) Bz averaged over the y direction, and e) x — Vx phase space. 
The zoomed in frame f) shows Bz in a 100 X 100 (c/ujpi)^ box at the shock front. 
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Fig. 15. — Panel a) shows the B^ magnetic field component for run §3. Panels b) and c) show lineouts of the By (blue) and B^ magnetic 
field components in the x and y directions along dashed lines. Frame d) shows the magnetic field for an identical run with a reduced 
transverse side (run Q4 with Ly = 4c/ujpi), and frame e) shows the lineout for By (blue) and Bz for the same run. 



run CJ3, which reproduces the physical parameters of the PIC run J^ in Figure 13 d (see Table 111). At times t = 40w^; 
(black lines) and t = 87 u!~^ (red line) the spectra from the two runs are very similar, and correspond to a downstream 
thermal spread of Eth ~ 0.45 E„j, {^ 45% of the incoming flow bulk energy). At this point in time, although there are 
non-thermal particles, the power law is very steep and the energy efficiency is very low. Measurements of the particle 
spectra at the same distance from the shock front at later times are shown in Figure |13^ , and shows the power-law 
part of the distribution growing in energy as time elapses. This is a typical signature of particle acce leration due 



to a diffusive process. The complete fit to the particle spectrum at t = 400 w^ is shown in Figure 13:, and yields 
Eth ~ 0. 43 E up, an injection energy of Einj '^ 1.45 E^p, and a power-law index of a ~ 1.9. 

Figure [14] shows the global structure for the Ma = 4.7 6 = 0° shock, simulated with the hybrid code, where a 
much larger L^ = 2000c/wpi Ly = lOOc/wpi box was used. It is evident in the zoomed in box for the B^ field 
component. Figure [T4ji, that there is a two-dimensional structure to the waves in the shock transition region. Also 
this wave structure is limited in size in the x direction, and defines the region where the downstream plasma is not 
yet completely thermalized. If we now consider a scenario where the transverse simulation size is of the same order 
of magnitude of the upstream ion Larmor radius r^i = 4.7c/ajpi, as it is i mpli citly done in ID simulations, again a 
coherent and extended transition region is formed. This is shown in Figure |15[ which compares the wave structures 
around the shock front for runs G3 (Ly = lOOc/wpi) and G4 (Ly = 4c/a;pi). 

These results show that numerical simulation parameters, and in particular the simulation box sizes, must be 
chosen with care. The transverse box size must be at least comparable to the typical wavelength observed for waves 
propagating in the x direction (parallel to the background magnetic field). To this requirement, we must add that the 
X box size must be such as to include a significant fraction of the pre-shock, including the highest energy particles. This 
requirement translates into very long box sizes for quasi-parallel (or strictly parallel) shocks, which greatly increases 
the cost of running these simulations for long times t > GOOoj^TJ . Thus, for higher Mach numbers (increasing V, while 



keeping B constant) larger x box sizes are needed, which justifies the L^ 
chosen for higher Mach numbers. 



5000c/cL'pi and L, 



8000 c/wpi box sizes 
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